#!/usr/bin/python

"""
Author: Jimmy Saw
Date: 07-03-2011
This program parses 454 paired-end reads aligned against an assembled genome.
Needs MUMMER output file generated like this:
nucmer genome.fasta 454reads.fasta -p file (Need to align reads binned as at least one having Gloeobacter)
show-coords -r -c -l -T file.delta > file.coords
Usage: dissertation_ParseOverlappingMatePairs.py file.coords > file.matepairs.txt
Usage: dissertation_ParseOverlappingMatePairs.py glbkl_vs_AtleastOneAndSingletons.coords > overlapping_clonepairs.txt

"""

import sys
import re

mummer_file = sys.argv[1]
mf = open(mummer_file, "rU")
mfl = mf.readlines()

left = re.compile('(\w+)_L')
right = re.compile('(\w+)_R')

paired = []
unpaired = []
allreads = []
reads_dict = {}
readnames = []

for line in mfl[4:]:
    strand = ""
    l = line.split('\t')
    ref_start = int(l[0])
    ref_stop = int(l[1])
    qry_start = int(l[2])
    qry_stop = int(l[3])
    if qry_stop > qry_start:
        strand = "+"
    else:
        strand = "-"
    per_id = float(l[6])
    qry_cov = float(l[10])
    qry_name = l[12].rstrip()
    readnames.append(qry_name)
    reads_dict[qry_name] = ((ref_start, ref_stop, qry_start, qry_stop, strand, per_id, qry_cov))

mate_pairs = []

for read in readnames:
    if left.match(read):
        pattern = left.match(read)
        mate_pair = pattern.group(1) + "_R"
        if read in reads_dict and mate_pair in reads_dict:
            if reads_dict[read][4] == "+":
                if reads_dict[mate_pair][4] == "-":
                    #calculate distance
                    if reads_dict[mate_pair][1] > reads_dict[read][1]:
                        pair_distance = reads_dict[mate_pair][1] - reads_dict[read][1]
                        ref_start_read = reads_dict[read][0]
                        ref_start_mate = reads_dict[mate_pair][1]
                        mate_pairs.append((ref_start_read, ref_start_mate, pattern.group(1), "L->R", pair_distance))
                    else:
                        pair_distance = reads_dict[read][1] - reads_dict[mate_pair][1]
                        ref_start_mate = reads_dict[mate_pair][0]
                        ref_start_read = reads_dict[read][1]
                        mate_pairs.append((ref_start_mate, ref_start_read, pattern.group(1), "R->L", pair_distance))
                else:
                    #bad alignment. shouldn't be in same orientation as read
                    pass
    elif right.match(read):
        pattern = right.match(read)
        mate_pair = pattern.group(1) + "_L"
        if pattern.group(1) in mate_pairs:
            pass #mate pairs already checked in previous statement
        else:
            if read in reads_dict and mate_pair in reads_dict:
                if reads_dict[read][4] == "+":
                    if reads_dict[mate_pair][4] == "-":
                        #calculate distance
                        if reads_dict[mate_pair][1] > reads_dict[read][1]:
                            pair_distance = reads_dict[mate_pair][1] - reads_dict[read][1]
                            ref_start_read = reads_dict[read][0]
                            ref_start_mate = reads_dict[mate_pair][1]
                            mate_pairs.append((ref_start_read, ref_start_mate, pattern.group(1), "R->L", pair_distance))
                        else:
                            pair_distance = reads_dict[read][1] - reads_dict[mate_pair][1]
                            ref_start_mate = reads_dict[mate_pair][0]
                            ref_start_read = reads_dict[read][1]
                            mate_pairs.append((ref_start_mate, ref_start_read, pattern.group(1), "L->R", pair_distance))
                    else:
                        #bad alignment. shouldn't be in same orientation as read
                        pass
    else:
        #singletons
        #unpaired.append(read)
        pass

mate_pairs.sort()

#Iterate over a given segment of the genome to sort mate pairs
tmpstr = mfl[5].split('\t')
genome_size = int(tmpstr[7])

#Print the mate pairs
for pair in mate_pairs:
    if 5000 < pair[4] < 12000: #only print pairs with expected insert size
        print pair[2] + "\t" + pair[3] + "\t" + str(pair[0]) + " -> " + str(pair[1]) + "\t" + str(pair[4])

mf.close()
